#!/usr/bin/perl -w																	#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use Getopt::Long;
use parse_bl2seq;
use db_parameters;
use File::Temp qw/ tempfile/;
use DBI;
use DBD::Pg;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

    align_sequences.pl tax1 tax2 file_ids
                          [--db_name rat_genome_new]
                          [--user postgres]
USAGE



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $dbname			= 'rat_genome_new';
my $username		= 'postgres';
my $userpassword	= '';
my $userhost		= 'fgu027';

GetOptions (

			'help'				=> \$help,
			'user=s'			=> \$username,
			'db_name=s'			=> \$dbname,
			);

die $usage if ($help);
die $usage unless @ARGV == 3;
my ($tax_name1) = shift @ARGV;
my ($tax_name2) = shift @ARGV;


#_________________________________________________________________________________________
#	sub declarations
#
sub do_align($$$$$);
#
#_________________________________________________________________________________________


#_________________________________________________________________________________________
#	Get DB parameters
#
use db_parameters;
my @db_param = get_db_parameters();
#
#
#_________________________________________________________________________________________




#_________________________________________________________________________________________

#     do_align

#        iterate through each ortholog pairs
#        get all proteins
#        blast2seq against each other and get smallest E-value
#        put results into

#_________________________________________________________________________________________


sub do_align($$$$$)
{
	my ($dbh, $gene_id1, $gene_id2, $tax_name1, $tax_name2) = @_;

	my @results = get_best_alignment($dbh, $gene_id1, $gene_id2, $tax_name1, $tax_name2,
						   {
							'seg'		=> 1,
							'coils'		=> 1,
							'e_value'	=> 1e-4,
							'blosom'	=> 80});
	return unless (@results);




	{
		mkdir ("yn00_tmp");
		open FH_SEQ, ">yn00_tmp/seq.tmp";
		print FH_SEQ "2 ".length($results[BL_CDNA_SEQ1]), "\n";
		print FH_SEQ "a\n", $results[BL_CDNA_SEQ1], "\n";
		print FH_SEQ "b\n", $results[BL_CDNA_SEQ2], "\n";
		close FH_SEQ;
		unlink "yn00_tmp/yn";
		open FH_YN00_CTL, ">yn00_tmp/yn00.ctl";
		print FH_YN00_CTL <<FILE;
seqfile = seq.tmp
outfile = yn
verbose = 0
noisy = 0
icode = 0
weighting = 1
commonf3x4 = 0
FILE
		close FH_YN00_CTL;
		chdir "yn00_tmp";
		system ("/share/tools/paml/yn00 < yn00.ctl > /dev/null");
		open FH_YN, "yn";
		chdir "..";
		while (<FH_YN>)
		{
			next unless /seq\.\s+seq./;
			<FH_YN>;
			$_ = <FH_YN>;
			my @values = split /[\s\+\-]+/;
			return unless @values == 12;
			return ($results[BL_BIT_SCORE],
					$results[BL_ALIGN_BEG],
					$results[BL_ALIGN_LEN],
					$results[BL_TOTAL_LEN],
					(int($results[BL_COVERAGE_] * 100.0 + 0.5)).'%' ,
					$results[BL_E___VALUE],
					$values[7],					# kaks
					$values[8],					# ka
					$values[3], 				# S
					$values[10]);				# Ks
		}
	}

}








#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#     main

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# connect to genome database
my $dbh = connect_to_panda();



print "gene_id1\tgene_id2\tbits\taln_beg\taln_len\ttot_len\t%cover\te-value\tKaKs\tKa\tS\tKs\tsteve_Ks\n";
#	my ($bitscore,
#		$beg,
#		$len,
#		$e_value,
#		$S
#		$ks) = do_align($gene_id1, $gene_id2, $tax_name1, $tax_name2);


for (<>)
{
	chomp;
	my ($gene_id1, $gene_id2, $old_ks) = split /\t/;
	if (!defined $old_ks)
	{
		print $_, "\n";
		next;
	}
	my (@results) = do_align($dbh, $gene_id1, $gene_id2, $tax_name1, $tax_name2);
	$results[9] = '' unless defined $results[9] and $results[9] < 50.0;
	$results[6] = '' unless defined $results[6] and $results[6] < 50.0;
	$results[7] = '' unless defined $results[7] and $results[7] < 50.0;
	@results = ('') x 13 unless (defined $results[0]);
#	my ($bitscore,
#		$beg,
#		$len,
#		$e_value,
#		$S
#		$ks) = do_align($gene_id1, $gene_id2, $tax_name1, $tax_name2);
	print join ("\t", $gene_id1, $gene_id2, @results, $old_ks), "\n";
#	last if $gene_id1 eq 'ENSG00000128645';
}

$dbh->disconnect();


